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Abstract 

In this paper, the outage probability and outage-based beam design for multiple-input multiple-output 
(MIMO) interference channels are considered. First, closed-form expressions for the outage probability 
in MIMO interference channels are derived under the assumption of Gaussian-distributed channel state 
information (CSI) error, and the asymptotic behavior of the outage probability as a function of several 
system parameters is examined by using the Chernoff bound. It is shown that the outage probability 
decreases exponentially with respect to the quality of CSI measured by the inverse of the mean square 
error of CSI. Second, based on the derived outage probability expressions, an iterative beam design 
algorithm for maximizing the sum outage rate is proposed. Numerical results show that the proposed 
beam design algorithm yields better sum outage rate performance than conventional algorithms such as 
interference alignment developed under the assumption of perfect CSI. 
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I. Introduction 

Due to their importance in current and future wireless communication systems, multiple-input multiple- 
output (MIMO) interference channels have gained much attention from the research community in recent 
years. Since Cadambe and Jafar showed that interference alignment (IA) achieved the maximum number 
of degrees of freedom in MIMO interference channels [2], there has been extensive research in devising 
good beam design algorithms for MIMO interference channels. Now, there are many available beam design 
algorithms for MIMO interference channels such as IA-based algorithms [3]-[5] and sum-rate targeted 
algorithms [3], [4], [6]-[9]. However, most of these algorithms assume perfect channel state information 
(CSI) at transmitters and receivers, whereas the assumption of perfect CSI is unrealistic in practical 
wireless communication systems since perfect CSI is unavailable in practical wireless communication 
systems due to channel estimation error, limited feedback or other limitations [10]. Thus, the CSI error 
should be incorporated into the beam design to yield better performance, and this is typically done under 
robust beam design frameworks. 

There are many robust beam design studies in the conventional single-user MIMO case and also in the 
multiple-input and single-output (MISO) multi-user case. In the MISO multi-user case, the problem is 
more tractable than in the MIMO multi-user case, and extensive research results are available on MISO 
broadcast and interference channels with imperfect CSI [1 1] — [13]; the outage rate region is defined 
for MISO interference channels in [11], and the optimal beam structure that achieves a Pareto-optimal 
point of the outage rate region is given in [12]. For more complicated MIMO interference channels, 
there are several pioneering works on robust beam design under CSI uncertainty [14]— [16]. In [14], the 
authors solved the problem based on a worst-case approach. In their work, the CSI error is modelled 
as a random variable under a Frobenius norm constraint, and a semi-definite relaxation method is used 
to obtain the beam vectors that maximize the minimum signal-to-interference-plus-noise ratio (SINR) 
over all users and all possible CSI error. In [15], on the other hand, the CSI error is modelled as an 
independent Gaussian random variable, and the beam is designed to minimize the mean square error 
(MSE) between the transmitted signal and the reconstructed signal at the receiver with given imperfect 
CSI at the transmitter (CSIT). 

In this paper, we consider the robust beam design in MIMO interference channels based on a different 
criterion. Here, we consider the rate outage due to channel uncertainty and the problem of sum rate 
maximization under an outage constraint in MIMO interference channels. This formulation is practically 
meaningful since an outage probability is assigned to each user and the supportable rate with the given 
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outage probability is maximized in practical systems. Here, we assume that the transmitters and receivers 
have imperfect CSI and the CSI error is circularly-symmetric complex Gaussian distributed. Under this 
assumption, we first derive closed-form expressions for the outage probability in MIMO interference 
channels for an arbitrarily given set of transmit and receive beamforming vectors, and then derive the 
asymptotic behavior of the outage probability as a function of several system parameters by using the 
Chernoff bound. It is shown that the outage probability decreases exponentially with respect to (w.r.t.) the 
quality of CSI measured by the inverse of the MSE of CSI, typically called the channel K factor [10] or 
interpreted as the Fisher information [17] in statistical estimation theory. In particular, it is shown that in 
the case of interference alignment, the outage probability can be made arbitrarily small by improving the 
CSI quality if the target rate is strictly less than the rate obtained by using the estimated as the nominal 
channel. Next, based on the derived outage probability expressions, we propose an iterative beam design 
algorithm for maximizing the weighted sum rate under the constraint that the outage probability for each 
user is less than a certain level. Numerical results show that the proposed beam design algorithm yields 
better sum outage rate performance than conventional beam design algorithms such as the 'max-SINR' 
algorithm [3] developed without the consideration of channel uncertainty. 

A. Related work 

The outage analysis for MIMO interference channels has been performed by several other researchers 
[16], [18]. In [16], the outage probability for a given rate tuple is computed under the assumption that 
the knowledge of the channel mean and covariance matrix are available, and transmit and receive beam 
vectors that minimize the power consumption for a given outage constraint are obtained. However, it is 
difficult to generalize this method of analysis to the case of multiple data streams per user, whereas our 
analysis includes the multiple data stream case. In [18], the outage probability and SINR distribution of 
each user in MIMO interference channels with the knowledge of channel distribution information are 
obtained under a particular transmit and receive beam structure of IA transmit beams and zero-forcing 
(ZF) receivers. On the other hand, our analysis can be applied to the case of general transmit and receive 
beam structures beyond IA and ZF. 

The probability distribution of a quadratic form of Gaussian random variables has been studied 
extensively in the statistics field [19]-[22] and in the communications area [23]-[25]. The most widely- 
used approach to obtain the probability distribution of a Gaussian quadratic form is the series fitting 
method [20], [21], [23], [26], which typically converges to the probability distribution of a Gaussian 
quadratic form from the lower tail first. However, the outage definition associated with robust beam design 
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for MIMO interference channels in this paper requires accurate computation of upper tail probabilities. 
The series expansion for the cumulative distribution function (CDF) obtained in this paper based on 
the integral form for the CDF in [25] and the residue theorem [22] is well suited to this purpose and 
converges to the upper tail first. Thus, the obtained series in this paper is more relevant for our outage 
analysis. For a detailed explanation of the derived series, please see Appendices B-C. 

B. Notation and organization 

We will make use of standard notational conventions. Vectors and matrices are written in boldface with 
matrices in capitals. All vectors are column vectors. For a matrix A, A H , \\A\\f and A(i,j) indicate 
the Hermitian transpose, the Frobenius norm and the element in row i and column j of A, respectively, 
and vec(A) and tr(A) denote the vector composed of the columns of A and the trace of A, respectively. 
For vector a, ||a|| and [a]j represent the 2-norm and the i-th element of a, respectively. I n stands for the 
identity matrix of size n (the subscript is included only when necessary), and diag(di, • • • , d n ) means 
a diagonal matrix with diagonal elements d\, ■ ■ ■ , d n . x ~ CJ\f(fi, S) means that the random vector x 
has the circularly-symmetric complex Gaussian distribution with mean vector and covariance matrix 
S. K = {1, 2, • • • , K}, l = yf—1, and \A\ denotes the cardinality of the set A. 

The paper is organized as follows. The system model and problem formulation are described in Section 
II. In Section III, closed-form expressions for the outage probability are derived, and the behavior of the 
outage probability as a function of several system parameters is examined by using the Chernoff bound. 
In Section IV, an outage-based beam design algorithm is proposed. Numerical results are provided in 
Section V, followed by the conclusion in Section VI. 

II. System Model and Problem Formulation 

In this paper, we consider a if -user time-invariant MIMO interference channel in which each transmitter 
equipped with N t antennas is paired with a receiver equipped with N r antennas, and interferes with all 
receivers other than the desired receiver. We assume that transmitter k transmits d (< mm(N t , N r )) 
independent data streams to receiver k paired with transmitter k. Then, the received signal at receiver k 
is given by 

K 

y k = H fcfc V fc s fc + HfcjVjS; + n fc , (1) 

i=l,i^fe 

where H^j is the N r x Nt channel matrix from transmitter i to receiver k; V, = [vfV-- ,vj d) ] is the 
N t x d transmit beamforming matrix with normalized column vectors at transmitter i, i.e., ||v- m ^|| = 1 
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for m = 1, ■ ■ ■ , d; and s, = [sf~\ • • • , s^] T is the tlx 1 symbol vector at transmitter i. We assume that 
the transmit symbol vector Sj is drawn from the zero-mean Gaussian distribution with unit variance, i.e., 
Sj ~ CA/"(0, I), and the additive noise vector is zero-mean Gaussian distributed with variance a 2 , i.e., 
rifc ~ CJ\f(0,a 2 T). We assume that the CSI available to the system is not perfect. That is, neither the 
transmitters nor the receivers have perfect CSI. For the imperfect CSI, we adopt the following model 

Hfej = Hki + Efej (2) 

for each (k, i) G JCxJC, where H^j is the unknown true channel, H^j is the channel state available to the 
transmitters and the receivers, and E^j is the error between the true and available channel information. 
For the CSI error E^j between the true and available channel information, we adopt the Kronecker error 
model which is widely used for MIMO systems to model the error correlation that may be caused by 
the transmit and receive antenna structure [10]. Under this model, the CSI error E^j is given by 

B ki = S^H^SJ 72 , with vec(H^ } ) ~CN(0,a 2 h I) for some a\ > 0, (3) 

where ~S t and S r are transmit and receive antenna correlation matrices, respectively, and the elements 
of are independent and identically distributed (i.i.d.) and are drawn from a circularly-symmetric 

zero-mean complex Gaussian distribution. The CSI uncertainty matrix E&j is a circularly-symmetric 1 
complex Gaussian random matrix with distribution vec(Efcj) ~ CAf(0, a 2 {Sf ® S P )) [10, p.90], and 
a\ is the parameter capturing the uncertainty level in CSI. We assume that the E^'s are independent 
across transmitter-receiver pairs (k,i). To specify the quality of CSI and signal reception, we define two 
parameters 

K {ki) _ = llHfcjjll _ ||Hfci||j, ( fe ) , = IjHfcfejjg 

Ch ' E{||E fei |||} ^ tr(5] T^ Sr) • a 2 ■ 

(ki) 

K ch is the channel K factor defined as the ratio of the power of the known channel part to that of the 
unknown channel part, representing the quality of CSI [10], and is the signal-to-noise ratio (SNR) 
at receiver k since V^, and are normalized in our formulation. Hereafter, we will use 'H to represent 
the collection of channel information {H-ki, ^t, 5] r } known to the transmitters and receivers. By using 
the receiver filter uj™^ || = 1), receiver k projects the received signal in (1) to recover the 

desired signal stream m: 

4 m) = ( u i m) ) H yfc = (ui m) ) H ( (H fcfc + E fcfe )V fe s fe + (Hfei + E fcl )V A + n k j . 

'The circular symmetry of a random matrix in form of AZB with constant matrices A and B and a circularly-symmetric 
complex Gaussian matrix Z can easily be shown by a similar technique to that used in the Appendix A. 
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We assume that the design of the transmit beamforming matrices {Vk,k G /C} and receive filters 
{Ufc = [u[ ,u^],k € /C} is based on the available CSI H. This model of beam design and 
signal transmission and reception captures many coherent linear beamforming MIMO schemes including 
interference alignment and sum rate maximizing beamforming schemes [3], [6], [27] in which transmit 
and receive beamforming matrices are designed based on available CSI at transmitters and receivers. 
Under this processing model, the SINR for stream m of user k is given by 

SINR^ m) |^= (4) 

|(ui m) ) g H fcfc vi m )| 2 

|(ui m V E fcfc v£ m > P + \("l m) ) H (H fefc + E fcfc )v« P + Ej =1 l(4 m y (H« + E fei )vf | 2 + ^ ' 

where the numerator of the right-hand side (RHS) in (4) is the desired signal power, and the first, 
second, third and fourth terms in the denominator of the RHS in (4) represent the interference purely 
by channel uncertainty, inter-stream interference, other user interference and thermal noise, respectively. 
(Here, the dependence of SINR on % is explicitly shown. Since the dependence is clear, the notation 
|^ will be omitted hereafter.) Because the {Efcj} are random, SINR^ is a random variable for given 
% and {Vfc(7-£), \5 k {H), k G /C}. Thus, an outage at stream m of user k occurs if the supportable rate 
determined by the received SINR (4) is below the target rate R k m \ and the outage probability is given 
by 

Pr{outage} = Pr |log 2 (l + SINR<f } ) < ijjf } . (5) 

By rearranging the terms in (4), the outage event can be expressed as 

K d , (m)|2 

E E 4T m 4T j) > ' U %^ V : ' - - 2 == r, (6) 

i=l j=l 2 k - 1 



where 



u^ m)H E fcfc v^ m) , i = k and j = m, 



• {m)H (H ki + B ki )^\ otherwise. 



Since the {E^} are circularly-symmetric complex Gaussian random matrices, {X^\i = 1, • • • ,K,j = 
1, • • • , d} are circularly-symmetric complex Gaussian random variables, and the left-hand side (LHS) of 
(6) is a quadratic form of non-central Gaussian random variables. To simplify notation, we will use 
vector form from here on. In vector form, (6) can be expressed as 

x (rn)H x (rn) > ^ (g) 
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where X 



(m) 



[y( ml ) 



y(m<!) y(ml) 



, xfe ] T . The elements of the mean vector /LtV 



( m )c_ 



E{x[ m) }) of xjf° are given by 



[Vk \(i-l)d+j 



0. 



k, j = m, 



u[ m)H H fci v u; , otherwise, 



(9) 



for i = 1, • • • , K and j = 1, • • • , d, and the covariance matrix El of x["^ is given by a block diagonal 
matrix, since {E^j, i = 1, • • • , if} are independent for different values of i, i.e., 



4 m) - E{(xl m) - E{X< m >})(X< m > - E{X< m >})*} = ding(E^, 
where the d x d sub-block matrix Sl? i s given by 



(10) 



Wife 



(m)If 



(l)ff T (d) 



(2)H„ (1) 
(2)if^ (2) 



(2)H_, (d) 



v?^E t v« 

2 1 % 

7. & 1. 



{d)H^ (d) 



(11) 



for each i. (The proof of (11) is given in Appendix A.) In the following sections, we will derive closed- 
form expressions for (5), investigate the behavior of the outage probability as a function of several 
parameters, and propose an outage-based beam design algorithm. 



III. The Computation of the Outage Probability 

In this section, we first derive a closed-form expression for the outage probability in the general case 
of the Kronecker CSI error model, and then consider special cases. After this, we examine the behavior 
of the outage probability as a function of several important system parameters based on the Chernoff 
bound. 



A. Closed-form expressions for the outage probability 

For a Gaussian random vector X ~ CA/"(/x, X) with the eigendecomposition of its covariance matrix 
£ = *A* H , the CDF of X H QX for some given Q is given by [25] 

i roo r(t,u)+f3) -c 

Pr{X H QX < T } = — - TWT^du (12) 

1 ^ - s 2ttJ_ 00 luj + P det(I + («j + /3)Q) 

for some /3 > such that I + /3Q is positive definite, where Q = A^Sf^Q* A 1/2 , x = A~ 1/2 *^/i 

and c = x H (j- + m+ff Q^ 1 ) ^ ^ mm nere on > we wn l derive closed-form series expressions for the 

CDF of the outage probability in several important cases by applying the residue theorem used in [22] 
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to the integral form (12) for the CDF. First, we consider the most general case of the Kronecker CSI 
error model. The outage probability in this case is given by the following theorem. 

Theorem 1: For given transmit and receive beamforming matrices {V^ = [v^ , • • • , vi ]} and {U^ = 
[uj^\ • • • , u[ ]} designed based on % = {Hfci, 5^, S r }, the outage probability for stream m of user k 
with the target rate under the CSI error model (2) and (3) is given by 

Pr {outage} = Pr{log 2 (l + SINR<; m) ) < R™} 

where r is given in (6); {Aj, i = 1, • • • , k} are all the distinct eigenvalues of the Kd x Kd covariance 

matrix Ejj™^ in (10) with eigendecomposition = *[ m ^A[ ^1 , Ki is the multiplicity 2 of the 

(i) 

eigenvalue A^; xl lS the element of vector 

(m) / * — i»T,( m )^ ( m ) /1 ^\ 

Xfc : = (H ) 2 *fc Mfe ( 14 ) 

corresponding to the j-th eigenvector of the eigenvalue Aj (1 < j < Ki), i.e., it is the j-th element of 
(AiI Ki ) _ 2^r^ ff ^ m \ (^"^ is a A'ti x «j matrix composed of the eigenvectors of sl"^ associated with 

ffi(s) = 7TT 7 T7^ ; (15) 

and gl n \s) is the n-th derivative of g^s) w.r.t. s. 

Proof: By using (12) and the facts Q = I and X^ m ^ ~ CA/"(/x^ m \ St ) in this case, we obtain the 

outage probability for stream m of user k in an integral form as 

, srr , \ 1 fP+too p sT „ _ E" = i Ix?'! 2 ) 

P r{ X««X« > t} = , - _ / _ • nt](w *, (.6. 



where s = (3 + to; (/? > 0). The outage probability (16) can be expressed as a contour integral: 

1 f e ST e~ ^'=' ' x 

2^ Jc "« nti( i +^) K 



Pr{X^xf) > r} = -^^T—^— ds, (17) 



=:F(5) 

where C is a contour of integration containing the imaginary axis and the whole left half plane of the 
complex plane. By the residue theorem, the sum of the residues at singular points of F(s) which do not 
have positive real parts yields the contour integral in (17) times 2ttl. It is easy to see that the singular 

2 Since S^ m ' > is a normal matrix, we have Kd = Ys7=i K,i - 
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points of F(s) are s = and s = — 1/Aj, i = 1, • • • , k. Since s[ ^ are all positive-definite, s£ is 
positive definite and Aj > for all i. So, the outage probability is given by 

K 

Pr-joutage} = 1 - (Res F(s) + ^ Res F(s)j . (18) 
It is also easy to see from (17) that the residue of F(s) at s = is Resi ? (s) = 1. To compute Res F(s), 

s=0 s=— 1/Xt 

for each i we introduce Gi(s) defined as 

/ 1 \ p r ( s ~ 1 / A P ^p=i i+i„(.-i/j i )^?=i «p I I 

Gi(a) := F ( s — -- 



KJ s-\l\ VC p=1 {\^\{s-\l\))^ 



I 1 /i 1 X i s -' V"! v 0) l 2 — V Apts-l/AJ y^Kp I ( 9 )|2 

3 T(S — 1/Ai) „ Z^j=l IXj I Z^p^i l + ApCs-l/Aj) ^9=1 IXp 



=:/i 

e -(^+E^lx? ) l») x e ^ E; ' llX ' 3) ' 2 x ( eTS x . 



-AW =*(») 

Now, the residue of at s = —1/Ai is transformed to that of Gj(s) at s = 0. The Laurent series 
expansion of fi(s) and the Taylor series expansion of gi(s) at s = are given respectively by 

/ 0')|2 > 



i ^ i , . , ^ i 

^5J^l A. iS ) and *w = E^ 

7 n— \ / n— 



(A**) 

By multiplying the two series and computing the coefficient of 1/s, we obtain the residue of Gi(s) at 
s = as 

' n=Ki — 1 \ / 

for each i Finally, substituting the residues into (18) yields (13). ■ 

To compute (13), we need to compute {A^}, {xP } an ^ tne higher order derivatives of gi(s). The first 
two terms are easy to compute since they are related with the mean vector of size Kd and the covariance 
matrix of size KdxKd. Furthermore, the higher order derivatives of gi(s) can also be computed efficiently 
based on recursion. (Please see Appendix C-A.) Note that in the case that the elements in (3) have 
difference variances, (13) is still valid since the difference variances only change the covariance matrix 
(10) and the outage expression depends on the covariance matrix (10) through {Aj} and {xP }• 

Next, we provide some useful corollaries to Theorem 1 regarding the outage probability in meaningful 
special cases. First, we consider the case in which a subset of channels are perfectly known at receiver k, 
i.e., Hfcj = Hfcj and E^j = for some i 6 1C. This corresponds to the case in which channel estimation 
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or CSI feedback for some links is easier than that for other links. For example, the desired link channel 
may be easier to estimate than others. The outage probability in this case is given by the following 
corollary. 

Corollary 1: When perfect CSI for some channel links including the desired link is available at receiver 
k, i.e., Hfcj = for i G C /C, the outage probability for stream m of user k is given by 

Pr{outage} = Pr{log 2 (l + SINR^ m) ) < R^} 

1=1 1 n=Ki— 1 \ / 

where r' is defined below; {Aj, i = 1, • • • , k'} is the set of all the distinct eigenvalues of the covariance 
matrix (10); Ki is the multiplicity of Aj, satisfying (K — | T^Dd = 2~2i=i K ^ Xi i s given in (14); and 



pynl-V (^-1/A0A P ^ Kp I (<?) |2 

LX P \ l^p^i l + f s -l/A,U„ 2^q=l \Xp 



5i,i(s) = -j- -, -■ (22) 

n p# , (l + («-l/A i )A p ) 



Proof: When CSI for some links including the desired link is perfect, the outage event at stream m 
of user k is given by 



|u<"^H fcfc vl m) | 2 ^ 



since = for i G T^. Thus, in this case the outage event is expressed in a quadratic form as 
follows: 

d I (m)Hxr ( m )|2 ^ 

EE4r )g 4r j) > |Ufc fl( ) fcfcVfc 1 - EEi«Wi 2 - £ lui^iw^-^w, 

»ex-j=i 2 k -1 igTfcJ - =1) ieXfc: 

(23) 

and we have = for all i G (See (7)). The size of xl m ^ now reduces to (K — |Tfc|)d, and 

the size of the covariance matrix E^ 1 i s ~~ |f fcl)^ x ~~ With the new threshold r', the 

same argument as that in Theorem 1 can be applied to yield the result. ■ 

Thus, when perfect CSI is available for some links, the order of the distribution is reduced under the same 
structure. Next, consider the specific beam design method of interference alignment and the corresponding 
outage probability, which can be obtained by Corollary 1 and is given in the following corollary. 

Corollary 2: When the desired channel link is perfectly known (i.e. k G T^) and {V^} and {U^} are 
designed under I A based on T-L, the outage probability for stream m of user k is given by 

Pr{outage} = -£-L e -£ 1 glf^O). (24) 

i=1 A i \ K i L >- 
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Proof: First, express the random term in (23) as J2ieri Ej=i X^ i)H X^ j) = (X^ m) )^X^ m) . 
When the beam is designed under IA based on %, we have E{Xjj. } = since u^^HjfcfVj = for 
all i e JC\{k} D T%, j = 1, ■ ■ ■ , d. (See (9).) Hence, xj; m) = and thus X ? = for all i and j. (See 
(14).) Then, the terms in the infinite series in (21) are zero for all n > m — 1 from the fact that 0° = 1 
and 0! = 1, and the result follows. ■ 



The outage probability for single stream communication is given in Corollary 3. 

Corollar 
is given by 



Corollary 3: When d = 1 and all eigenvalues of are distinct, the outage probability for user k 



Pr{outage} = Pr{log 2 (l + SINR fc ) < R k } 



/A -£(- n 

n=0 ^ 

XpCa-l/Xj) 



n— v 

A p (s-1/A,) 



i U"»(0), 



where gi(s) in (15) reduces to gi(s) = —ztt\ ' ~ 7 — — "' — \~ • (Here, we have omitted the stream 

superscripts since the stream index is unique.) 

Proof: Since all eigenvalues are assumed to be distinct, there are k = K eigenvalues with m = 1 
for all i. Substituting these into Theorem 1 yields the result. ■ 



Now, let us consider a simpler case for d = 1 with no antenna correlation. In this case, the outage 



probability is given as an explicit function of the channel uncertainty level a?, and it is given by the 



following corollary to Theorem 1. 

Corollary 4: When d = 1 and there is no antenna correlation, the outage probability is given by 

where Xk = E{X fc }/a h and g(s) = g _ e 1 " CT . ■ 

Proof: In this case, an outage at user k occurs if and only if X^X^ > ^ Ufc 2 ^ fc fc ^' c a 2 . Now, the 

covariance matrix of X^ is p^Ik (see (10) and (11)), and thus there is only one eigenvalue a\ with 
multiplicity K. Moreover, Xk = ^{^-k}/^h from (14) since = I and = a%I. By substituting 
these into Theorem 1, the outage probability (26) is obtained. ■ 

B. The behavior analysis of the outage probability based on the Chernoff bound 

The obtained exact expressions for the outage probability in the previous subsection can easily be 
computed numerically, and will be used for the robust beam design based on the outage probability 
in Section IV. Before we address the outage-based robust beam design problem, let us investigate 
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the behavior of the outage probability as a function of several parameters. Suppose that transmit and 
receive beam vectors {vj|, m \u^} are designed by some known method based on %. For the given 
beam vectors, as seen in the obtained expressions, the outage probability is a function of other system 
parameters such as the known channel mean {H^j}, the noise variance a 2 , the channel uncertainty level 
a 2 ., the antenna correlation E t and S r , and the target rate fl£ . Here, the dependence on Hjy-, a 2 and 
R k m ^ is via the threshold r(Hfcfc, a 2 , R^), and the dependence on a\, S t , S r and {Hki,i ^ k} is 
via *J; m) (4 m) K, St. S r ), E{X^ m) }(H fci )) and the eigenvalues of sJ^K, S t , S r ). This complicated 
dependence structure makes it difficult to analyze the properties of the outage probability as a function 
of the system parameters. Thus, in this subsection we apply the Chernoff bounding technique [17] to the 
tractable 3 case of d = 1 to obtain insights into the outage probability as a function of several important 
parameters. When d = 1, the outage event is expressed as 

Pr{xf X fc > r = l ^t k _\f " - 2 } = Pr { * ^ < 27 > 

Since Ej-i,--- , Efcx are independent and circularly-symmetric complex Gaussian random matrices, 
Xki, • • • , XkK are independent and circularly-symmetric complex Gaussian random variables. (See (7).) 
Thus, the term on the LHS in the second bracket in (27) is a sum of independent random variables, and 
the Chernoff bound can be applied to yield 

Pr{XfX fe >r}< e- TS Uli^{e slX ^} (28) 

for any s > 0. The moment generating function (m.g.f.) of \Xki\ 2 (Xki ~ a "ki)) i s given by 

E {yl* M | 2 } = exp f^l) for 8 < 1/al, where m = 0, fx hi = uf H fei Vi for i + k, and a\ 



a^(u^S r .Ufc)(v^S i Vj). (See (7,9,11).) Therefore, the Chernoff bound on the outage probability is given 
by 

A 



Pr{X^X fc >r}< e— J] T^" ex P fr^) 




K , 



(29) 



for < s < minj {l/o-|J. Now, (29) provid es a tool to analyze the behavior of the outage probability 
as a function of several important parameters. The most desired property is the behavior of the outage 

3 In certain cases of d > 1, Chernoff bound can still be obtained when each element in Xj_ m ' is independent of the others. 
Such cases include the case that there is no antenna correlation and the transmit beam vectors are orthogonal as in the IA beam 
case. In this case, similar results to the case of d = 1 are obtained. 
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probability as a function of the channel uncertainty level. This behavior is explained in the following 
theorem. 



Theorem 2: When d 
given by 



1, as g\ — > 0, the outage probability decreases to zero, and the decay rate is 



Prjoutage} < e Cl • exp(— c^jo^) 



(30) 



for some c\ and C2 > not depending on a\, if the target rate and the designed transmit and receive 
beam vectors {vfc,Ufc} satisfy 

/ \ 



Ru < Ri. 



log 2 



1 + 



\ 



ife: 



mi 



(31) 



/ 



tr(E r )t7-(E t ) 

Proof: (29) is valid for any s G (0, min^l/cr^}). So, let s = l/oj[tr(E t )tr(E r ) (< min^l/cr^} 
since ||v fc || = ||u fc || = 1 and = fj^(u^E r .u fc )(vf^E t v i ) < cr^tr(Ei)tr(E r ) for all i). Then, the 
exponent in (29) is given by 



K 



^tr(E t )tr(E r ) 



K 

E 



i - 



(ufS r u fc )(vfS tVi ) 



tr(S t )tr(S r ) 
l/^| 2 



a£ \tr(E t )tr(E r ) ^ (ujf S r u fc )(vf E tVi ) - tr(E t )tr(E r ; 



if 

E ^(ufl S rUfc )(vf E tVi ) - a2tr(S t )tr(S r ) 



(ufS r .u fc )(vfS t v t ) 
tr(E t )tr(E P ) 



>-f>g 



1 - 



(=:<*) 

.-ffTT. . „. 12 //-Ofl* 



Now, substituting r = |u^HfcfcVfc| 2 /(2 Hk — 1) — a into the inequality c 2 > yields (31). ■ 

Theorem 2 states that the outage probability decays to zero as the CSI quality improves, more precisely, 
it decays exponentially w.r.t. the inverse of channel estimation MSE (or equivalently w.r.t. the channel 
K factor), if the target rate is below R^. In the Fisherian inference framework, the inverse of estimation 
MSE is information. Thus, another way we can view the above is that the outage probability decays 
exponentially as the Fisher information for channel state increases, if the target rate is below a certain 
value. So, the outage probability due to channel uncertainty is another case in which information is 
the error exponent as in many other inference problems. In certain cases, the condition (31) can be 
simplified considerably. For example, when interference-aligning beam vectors based on % are used at 
the transmitters and receivers, we have fi^i = u^H^Vj = for i / k in addition to /j,^ = 0, and the 

/ lu^H v I 2 \ 

condition is simplified to R^ < log 2 ( 1 + k ^ I . Thus, in the case of interference alignment the 
outage probability can be made arbitrarily small by improving the CSI quality if the target rate is strictly 
less than the rate obtained by using as the nominal channel. Next, consider the outage behavior as 
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the effective SNR, T e ff := |uj^HfcfcVfc| 2 /<r 2 , increases. Since the two terms determining the effective 
SNR are contained only in r, it is straightforward to see from (29) that 

Pr{outage} < c 3 exp (-c 4 r e// ) , (32) 

for some C3 and C4 = sa 2 /(2 Rk — 1) > not depending on T e ff- Finally, consider the case in which 
the target rate decreases. One can expect that the outage probability decays to zero if the target rate 
decreases to zero. The decaying behavior in this case is given in the following theorem. 

Theorem 3: When d = 1, as — > 0, the outage probability decreases to zero, and the decay rate is 
given by 

Pr{outage} < c 6 exp (-p^y) = c 6 exp (- ^^J (33) 

for some C7, c' 7 > not depending on R^. The last equality is when R^ is near zero. 

Proof: Let s be any positive constant contained in an interval (0, 1/ maxj{cr 2 1 (u^X! r u/ c )(v^X t Vj)}). 
Then, the exponent in (29) becomes 

i=1 i=1 S<7 h\ U k ■^rU fc J(V i 2j t Vi) - 1 

\ ✓ 

(=:ob) 

/|ufH fcfc v fe | 2 _ 2 \_ |ufH fcfc v fc | 2 _ 



2^-1 5 " c s-- a^.i 



Hence, the Chernoff bound is given by Prjoutage} < cq exp ( — s l u fc^ Vfc j — Cg eX p |_ c j j f or 



2^-1 J ~ ^ ^ \-R k +o(R k ) / 

some c'j > 0. The last equality is when R^ is near zero. In this case, we have 2 Rk — 1 = (log 2)Rk + o(Rk) 
by Taylor's expansion. ■ 

IV. Outage-Based Robust Beam Design 

In this section, we propose an outage-based beam design algorithm based on the closed-form expres- 
sions for the outage probability derived in the previous section. Our assumption is that % is given for the 
beam design, as mentioned earlier. Suppose that transmit and receive beamforming matrices {Vfc,Ufc} 
are designed by using any available beam design method based on %. Based on the designed {V^, U^} 
and known {H,a 2 }, one can compute and use a nominal rate for transmission. Since % is not perfect, 
however, an outage may occur depending on the CSI error if the nominal rate is used for transmission. Of 
course, the outage probability can be made small by making the transmission rate low or by improving 
the CSI quality, as seen in Section III-B. However, these methods are inefficient sometimes since we 
may have limitations in the CSI quality or need as high rate as possible for given %. Further, in many 
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wireless systems the target outage probability for transmission is determined and the data transmission 
is performed under such an outage constraint. Thus, we here consider the beam design problem when 
the outage probability is given as a system parameter. In particular, we consider the following per-stream 
based beam design problem to maximize the sum e-outage rate for given T-L: 



maximize 

{vi m) },{u<™>} 



E £ (34) 

k=l m=l 

subject to Pr{log 2 (l + SINRjf 1^) < R^} < e (35) 
Il4 m) || = ||v£ m) || = 1, VfceJC, m = l,-.- ,d, (36) 

where the e-outage rate for stream m of user k is the maximum rate satisfying (35). Like other beam 
design problems in MIMO interference channels, the simultaneous joint optimal design for all transmit 
and receive beam vectors for this problem also seems difficult. Hence, we propose an iterative approach to 
the above sum e-outage rate maximization problem. The proposed method is explained as follows. In the 
first step, we initialize {v[ } and {u[ } properly (here a known beam design algorithm for the MIMO 
interference channel can be used), and then find optimal rate-tuple (R^\ • • • , R^ , R^\ " * * > Rjc) tnat 
maximizes the sum for given {vjj. m \ ujj."^} under the outage constraint. This step is performed based on 
the derived outage probability expressions in the previous section. Since designing each R^ does not 
affect others, this step can be done separately for each R^ ■ Since the outage probability for stream m 
of user k increases monotonically w.r.t. R^\ the optimal R^ in this step is the rate with the outage 
probability e. In the second step, for the obtained rate-tuple and receive beam vectors {u[ m ^} in the first 
step, we update the transmit beam vectors {v^ m ^} to minimize the maximum of the outage probabilities 
of all streams and all users. (Since the outage probabilities of all streams of all users are e at the end 
of the first step, this means that the outage probability decreases for all streams and all users.) Here, we 
apply the alternating minimization technique [28] to circumvent the difficulty in the joint transmit beam 
design. (The change in one transmit beam vector affects the outage probabilities of other users.) That 
is, we optimize one transmit beam vector while fixing all the others at a time. We iterate this procedure 
from the first stream of transmitter 1 to the last stream of user K until this step converge. In the third 
step, we design the receive beam vector vS™^ to minimize the outage probability at stream m of user k 
with the rate-tuple determined in the first step and {v^} determined in the second step for each (k, m). 
This optimization can also be performed separately for each stream of each user since the receiver filter 
for one stream does not affect the performance of other streams. Finally, we go back to the first step 
with the updated transmit and receive beam vectors (in the revisited first step, the rate for each stream 
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will be increased by increasing the outage probability upto to e again), and iterate the procedure until the 
sum e-outage rate does not change. We have summarized the sum outage rate maximizing beam design 
algorithm in Table I. 

The Proposed Algorithm 



Input: channel state estimate H and allowed outage probability e. 

(m)i , r In 

r l } and K 



0. Initialize {v^ m ' ) } and {uj^ } as sets of unit-norm vectors properly. 



1. For given {Vfc} and {Ufc}, find (R^\ • • ■ ,Rr - ) that maximizes J2k=i El=i Rk while 
the outage constraint is satisfied. 

2. Update {V fe = [v^, ■ • • ,vjj. d) ]} for {R^ n) } and {u[ m) } given from step 1. 

•For pair fix {vjj. , k = 1, • • • , K, m = 1, • • • , d}\{v^'} and {Ufc} and solve 

= argminmaxPr{outage[ m ' ) }. (37) 

vgC N t k,m 

(Here, a commercial tool such as the matlab fminimax function can be used to solve (37) 
together with the derived outage expression.) 

•Iterate the above step from the first stream of transmitter 1 to the last stream of transmitter 
K until {Vi, ■ • • , Vr-} converges. 

3. For receiver 1 to K, obtain the receive filter uj: m) that minimize the outage probability of 
stream m of receiver k for given {V^} from step 2 and given R ™ from step 1. (Here, 
again a commercial tool such as the matlab fmincon function can be used together with the 
derived outage expression.) 

4. Go to step 1 and repeat the whole procedure until the algorithm converges. 

TABLE I 

THE PROPOSED ALGORITHM FOR SUM e-OUTAGE rate maximization with channel uncertainty 



Theorem 4: The proposed beam design algorithm converges. 

Proof: It is straightforward to see that the sum e-outage rate increases monotonically for each iteration 
of the three steps of the proposed algorithm. Also, the maximum sum rate is bounded by the rate with 
perfect CSI. Hence, the algorithm converges by the monotone convergence theorem for real sequences. 

■ 

V. Numerical Results 

In this section, we provide some numerical results to validate our series derivation, to examine the 
outage probability as a function of several system parameters and to evaluate the performance of the 
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Fig. 1. Comparison of two series expressions for the CDF of quadratic form of Gaussian random variables. X ~ 
CA/"([0.5,0.5,0.5, 0.5] T ,0.3I 4 ), Q = [1, 0.5, 0, 0; 0.5, 1,0,0; 0,0, 1,0; 0,0,0, 1], and j3 = 2 for Laguerre series expansion. 



proposed beam design algorithm. For given Sj, S r , and T^ k \ we first generated {Hfcj} randomly 

according to zero-mean Gaussian distribution, and then scaled H^j to yield ||H^|||, = N t N r for all 
(k, i). In this way, the channel K factor and the SNR were simply controlled by a\ and a 2 , respectively. 
After {Hfcj} were generated as such, we generated {E^j} according to (3) and the true channel was 
determined by (2) if necessary 4 . For simplicity, we used = K c h for all (k, i) and r( fc ) = T for all 

k. 

First, Fig. 1 compares the convergence behavior of the derived series in this paper with that of the series 
fitting method [20], [21], [23], [26] based on the Laguerre basis functions for a given set of parameters 
shown in the label of the figure. It is seen that indeed our series converges from the upper tail first 
whereas the series fitting method converges from the lower tail first. (For a proof of this in the identity 
covariance matrix case, please refer to Appendix C-B.) Note that the series fitting method yields large 
error at the upper tail distribution even with a reasonably large number of terms. With this verification, 
next consider the outage behavior as a function of several system parameters. 

4 The computation of the closed-form outage probability requires only the channel statistics and {Tiki} regarding the channel 
information, but for Monte Carlo runs we need to generate {E^}. 
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Target Rate [bits/channel use] 

Fig. 2. Outage probability versus the target rate Rk (K = 3, Nt = N r = 2d = 2, Et = E r = I, T = 15 dB. Transmit and 
receive beam vectors are obtained by the IIA algorithm in [3].) 

Fig. 2 shows the outage probability w.r.t. the target rate R^ for a given set {Hfc«} (randomly generated 
as above) with several different channel K factors, when K = 3, N t = N r = 2d = 2, H t = S r = I, 
T = 15 dB and the transmit and receive beam vectors were designed by the iterative interference alignment 
(IIA) algorithm [3]. The solid and dotted lines represent the result of our analysis, and the markers + 
and x indicate the result of Monte Carlo runs for the outage probability. The theoretical outage curves in 
Fig. 2 were obtained by using (21) with the first 38 terms in the infinite series. It is seen that our analysis 
matches the result of Monte Carlo runs very well. The dashed line shows the outage performance when 
K c h = oo, i.e., all transmitters and receivers have perfect CSI. In the case of K c h = oo, we have a sharp 
transition behavior across Ru m it determined by the SINR (4) with E^j = for all (k,i). It is seen that 
the outage performance deteriorates from the ideal step curve of K c h = oo, as the CSI quality degrades. 
The solid lines correspond to the outage performance for the finite values of K c h, when the CSI for all 
channel links is imperfect. It is seen that = 100 (20 dB) yields reasonable outage performance 
compared with the perfect CSI case in this setup. Note that the gain in the outage probability by knowing 
the desired link perfectly is not negligible. (See the dotted lines.) Fig. 3 show the outage probability w.r.t. 
the target rate for a given set {H^.j} with several different K c h, when K = 3, Nt = N r = 2d = 4, 
St = S r = I, T = 25 dB and the transmit and receive beam vectors were designed by the IIA algorithm. 
Similar behavior is seen as in the single stream case, i.e., the outage performance generally deteriorates 
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as K c h decreases. However, it is interesting to observe in the multiple stream case that sufficiently good 
but not perfect CSI quality yields better outage performance than the perfect CSI in the high outage 
probability regime. (See Fig. 3 (b).) This implies that in the multiple stream case the second term (i.e., 
the self inter-stream interference term) in the denominator of the SINR formula (4) is made smaller by 
Efcfc's being negatively aligned with than in the case of E^, = 0. However, this is not useful in 
system operation since the system is operated in the low outage probability regime. All the theoretical 
curves in Figures 3 (a) and (b) were obtained by (21) with the first 45 terms in the infinite series. Fig. 4 
shows the outage probability curves when the transmit and receive beamforming vectors are respectively 
chosen as the right and left singular vectors corresponding to the largest singular value of the desired 
channel and the other parameters are identical to the case in Fig. 2. A similar outage probability behavior 
to the previous case is observed. 




1 2 3 4 5 6 7 

Target Rate [bits/channel use] 



(a) CDF 



o 

it O.i 



0.2 
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« = 3 x 10 4 






+ K ch = 10'' 






• K ch = 10" 
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Target Rate [bits/channel use] 
(b) Residual error 



Fig. 3. Outage probability versus the target rate R k (K = 3, N t = N r = 2d = 4, S t = S r = I, T = 25 dB. Transmit and 
receive beam vectors are designed by the IIA algorithm in [3].) 



Next, the outage probability w.r.t. the channel K factor for a given set {H} for several values of the 
target rate Rf. is shown in Fig. 5, where the outage probability along the y-axis is drawn in log scale. 
(The same setup as for Fig. 2 was used and the IIA algorithm is used for the transmit and receive beam 
design. Here, (21) with the first 38 terms in the infinite series was used to compute the analytic curves.) 
As predicted by Theorem 2, the outage probability indeed decays exponentially w.r.t. the channel K 
factor (equivalently, w.r.t. the inverse of af). The exponent depends on the target rate R^, the higher the 
target rate is, the smaller the exponent is. This decaying behavior is also predicted in Theorem 2; the 
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exponent C2 in (30) is proportional to r, and r is inversely proportional to the target rate Rk- It is seen 
that the outage probability does not decay as K c h increases, if Rk is larger than Rumit- In addition to 
the exact outage probability, the Chernoff bound in this case is shown in Fig. 5 as the lines with dots 
and dashes. It is seen that the Chernoff bound is not very tight but the decaying slope is the same as 
that of the exact outage probability. 




Target Rate [bits/channel use] 

Fig. 4. Outage probability versus the target rate R k (K = 3, N t = N T = 2d = 2, S t = S r = I, T = 15 dB. Transmit and 
receive beam vectors are respectively chosen as the right and left singular vectors corresponding to the largest singular value of 
the desired channel matrix.) 

Figures 6 and 7 show the impact of antenna correlation on the outage probability. We adopted the 
exponential antenna correlation profile considered in [29], [30]. Under this model, the (i, j)-th element 
of the antenna correlation matrix Si (or S r ) in (3) is given by p\ % ~i\ where p G [0, 1] is a parameter 
determining the correlation strength. Since tr(S t ) = N t and tr(S r ) = N r for this exponential antenna 
correlation model, we have the same transmit and receive powers as in the case of no antenna correlation, 
i.e., = I and S r = I. Since the outage probability depends on {Hfcj} as well as on St and S r , 
we generated one hundred {H^} randomly in the way that we explained already, and averaged the 
corresponding 100 outage probabilities to see the impact of the error correlation only. Other aspects of 
the system configuration were the same as those for Figures 2 and 5. It is seen that the error correlation 
decreases the outage probability especially when the CSI quality is very bad, but the gain becomes 
negligible when the CSI quality is good. 
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Fig. 5. Outage probability versus K c h (K — 3, N t = N r = 2d = 2, S t = S r = I, T = 15 dB. Transmit and receive beam 
vectors are designed by the IIA algorithm in [3].) 
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Fig. 6. Average outage probability versus T (K = 3, JVt = N r — 2d — 2. Transmit and receive beam vectors designed by the 
IIA algorithm in [3].) 
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Fig. 7. Average outage probability versus T (K — 3, N t — N r — 2d — 2, Rk — 1.2. Transmit and receive beam vectors 
designed by the IIA algorithm in [3].) 



Finally, the performance of the proposed beam design algorithm maximizing the sum e-outage rate 
was evaluated. As reference, we adopted the max-SINR algorithm and IIA algorithm in [3]. Although 
the max-SINR and IIA algorithms were originally proposed to design beam vectors with perfect channel 
information, we applied the algorithms to design beam vectors by treating the imperfect channel H as 
the true channel. The e-outage rate of the max-SINR algorithm (or the IIA algorithm) is defined as the 
maximum rate that can be achieved under the outage constraint of e using the beam vectors designed 
by the max-SINR algorithm (or the IIA algorithm). Once {V&} and {Ufc} are designed by any design 
method for given Ht, S r and {Hfej}, the outage probability corresponding to the designed beam vectors 
is easily computed as a function of the target rate Rk from Theorem 1. Thus, for the beam vectors 
designed by the max-SINR and IIA algorithms as well as for those designed by the proposed design 
algorithm in Section IV, the e-outage rate Rk can easily be obtained. Figures 8 and 9 show the sum e- 
outage rate of the proposed beam design method averaged over thirty different sets of {Hfcj} for e = 0.1 
and e = 0.2, respectively, when K = 3, N t = N r = 2d = 2 and E t = S r = I for different K^s. (The 
outage probability expression (26) with the first 40 terms was used to compute the outage probability.) 
It is seen that the proposed algorithm outperforms the IIA and max-SINR algorithms in all SNR, and 
the max-SINR algorithm shows good performance almost comparable to the proposed algorithm at low 
SNR. However, as SNR increases, the performance of the max-SINR algorithm degrades to that of the 
IIA algorithm (the two algorithm themselves converge as SNR increases) and there is a considerable gain 
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Fig. 8. Sum e-outage rate for e = 0.1 (K = 3, N t = N r = 2d = 2, S r = E r = I) 
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Fig. 9. Sum e-outage rate for e = 0.2 (K = 3, N t = N r = 2d = 2, S t = S r = I) 
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VI. Conclusion 

In this paper, we have considered the outate probability and the outage-based beam design for MIMO 
interference channels. We have derived closed-form expressions for the outage probability in MIMO 
interference channels under the assumption of Gaussian-distributed CSI error, and have derived the 
asymptotic behavior of the outage probability as a function of several system parameters based on the 
Chernoff bound. We have shown that the outage probability decreases exponentially w.r.t. the channel K 
factor defined as the ratio of the power of the known channel part and that of the unknown channel part. 
We have also provided an iterative beam design algorithm for maximizing the sum outage rate based 
on the derived outage probability expressions. Numerical results show that the proposed beam design 
method significantly outperforms conventional methods assuming perfect CSI in the sum outage rate 
performance. 

Appendix A 
Proof of (11) 

The (p, q)-th element of El"? is given by 

E{(J!r&*> -E{X^})(X^ -E[xW>}) B } 

( = } E{(vf T ® ui m ^)vec(E fel )vec(E fc4 )^(v^ T ® uj^)*} 

( = } 4i^ )T ® uf^XE? ® S r )(vf ® nt )H f 

( = } oftv^S? « u£»>*E r )(v?>- ® ui™>), where vp = (vfT 

Here, (a) is obtained by applying vec(ABC) = (C T ® A)vec(B) to each of the two terms in the 
expectation, (b) is by E{vec(E fci )vec(E fci ) w } = cr^(Sf ® S r ), (c) and (d) are by (A <g> B)(C ® D) = 
(AC BD), and finally (e) is because v^ T Sf * and u^ m ^5] r u^ m ' > are scalars. ■ 

Appendix B 

Distribution of a Non-Central Gaussian Quadratic Form 

The contents in Appendices B and C are from the technical report WISRL-2012-APR-1, KAIST, "A 
Study on the Series Expansion of Gaussian Quadratic Forms". 
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A. Previous work and literature survey 

There exist extensive literature about the probability distribution and statistical properties of a quadratic 
form of non-central (complex) Gaussian random variables in the communications area and the probability 
and statistics community. Through a literature survey, we found that the main technique to compute the 
distribution of a central (or a non-central) Gaussian quadratic form is based on series fitting, which was 
concretely unified and developed by S. Kotz [20], [21], and most of other works are its variants, e.g., 
[23]. First, we briefly explain this series fitting method here. 

Consider a Gaussian quadratic form x^Qx, where x ~ CjV(/Lt, S) with size n and Q = Q H . The 
first step of the series fitting method is to convert the non-central Gaussian quadratic form into a linear 
combination of chi-square random variables: 

n n 

X^Qx = M z i + = Yl X ii R < Z i + 6 i) 2 + lm ( Z i + S i)% < 38 ) 
i=l i=l 

where z% ~ CAf(0, 2) for % = 1, • • • , n, and {Si, X{} are constants determined by Q, /x and E. 
Note that Re(zj) ~ jV(0, 1) and Re(zj) ~ A/"(0, 1). Thus, the non-central Gaussian quadratic form is 
equivalent to a weighted sum of non-central Chi-square random variables of which moment generating 
function (MGF) is known. The MGF of a weighted sum of n independent non-central \ 2 random variables 
with degrees of freedom 2rrii and non-centrality parameter /i? is given by 

-i n i n 2 n i 

•m=«p{- E*?+ ET^Mladb-. (39) 

i=l i=l % i=\ y 1 ' 

Note here that <&(— s) is nothing but the Laplace transform of the probability density function (PDF) of 
x^Qx or equivalently J27=i x i\ z i + ^| 2 - Now, the series fitting method expresses the PDF as an infinite 
series composed of a set of known basis functions and tries to find the linear combination coefficients 
so that the Laplace transform of this series is the same as the known $(— s). Specifically, let the PDF be 

oo 

#„(Q,/x,£;y) = ^c k h k {y), (40) 

fc=0 

where {h k (y) : k = 0, 1, • • • } is the set of known basis functions and {c k , k = 0, 1, • • • } is the set of 
linear combination coefficients to be determined. Here, to make the problem tractable, in most cases, 
the following conditions are imposed. First, the sequence {/ifc(y)} of basis functions is chosen among 
measurable complex- valued functions on [0, oo] such that 

oo 

|cfc||/ifc(g/)| < Ae by , y £ [0,oo] almost everywhere, (41) 

fc=0 
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Gaussian quadratic form 



Linear combination of 
Chi-square random variables 



Moment generating function 
(or Laplace transform) 



Equivalent series expansion 

Kotz (67) 
Nabar (05) 



Direct inversion 
Raphaeli (96) 



Fig. 10. Computation of the distribution of a Gaussian quadratic form 



where A and b are real constants. Second, the Laplace transform hk{s) of hk(y) has a special form: 



(42) 



where is a non-vanishing, analytic function for Re(s) > b, and rj(s) is analytic for Re(s) > b and 
has an inverse function. The first condition is for the existence of Laplace transform and the second 
condition is to make the problem tractable. Finally, with the pre-determined {hk(y)} with the conditions, 
the coefficients {c^} are computed so that 



£(Sn(Q, S; y)) = c kh(s) = *(-«)> 



(43) 



fc=0 



where £(•) denote the Laplace transform of a function. 

Widely used {/tfc(y)} for the series expansion of the PDF of a quadratic form of non-central Gaussian 
random variables is as follows [20], [21]. 

1. (Power series): h k (y) = (-l) k ■ 

2. (Laguerre polynomials): 



h k (y) = a^y/m- ^n/^l Q ^k'^W 2 ®' 



(44) 
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where g(n; y) is the central x 2 density with n degrees of freedom and L^ 1 ^ 2 (x) is the generalized 
Laguerre polynomial defined by Rodriges' formula 

K k\ dx k 

for a > 1 and a positive control parameter /3. 

For the detail computation of {q,.}, please refer to [20], [21], [26]. The whole procedure is summarized 
in Fig. 10. 
Reference group 1 
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Gaussian random vector in Laguerre polynomial. 
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Indian Journal of Statistics vol. 17, pp. 37 - 50, Jan. 1956. =>■ CDF for the indefinite quadratic form 
of central random variable. 
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[Imhof-61] J. P. Imhof, "Computing the distribution of quadratic forms in normal variables," Biometrika 
vol. 48, pp. 419 - 426, Dec. 1961. =4> Provides a numerical method of computing the distribution 
[Rice-80] S. O. Rice, "Distribution of quadratic forms in normal variables - Evaluation by numerical 
integration," SIAM J. Scient. Statist. Comput., vol. 1, no. 4, pp. 438 -448, 1980. =>■ Another numerical 
method of computing distribution. 

[Biyari-93] K. H. Biyari and W. C. Lindsey, "Statistical distribution of Hermitian quadratic forms 
in complex Gaussian variables," IEEE Trans. Inform. Theory, vol. 39, pp. 1076 - 1082, Mar. 1993. 
=> Series expansion of multi-variate complex Gaussian random variables. This paper deals with the 
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[Raphaeli-96] D. Raphaeli, "Distribution of noncentral indefinite quadratic forms in complex normal 
variables," IEEE Trans. Inf. Theory, vol. 42, pp. 1002 - 1007, May 1996. 
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B. The difference of our work from the previous works 

First, let us remind our outage event in MIMO interference channels. From equations (5), (6) and (7), 
we have 

{K d i (m)_ff-pr (m),2 'l 

E E 4T i)H 4T j) > |Ufc 2fl " fc * V i - a2 = :T \ ' (45) 

where xj?™^ is a non zero-mean Gaussian random variable. Note that the outage probability is an upper 
tail probability of the distribution of the Gaussian quadratic form Ylj=i X^^ H However, as seen 
in Fig. 1 1 , the most widely-used series fitting method explained in the previous subsection yields a good 
approximation of the distribution at the lower tail not at the upper tail. The discrepancy between the 
series and the true PDF is large at the upper 5 tail for a truncated series. On the other hand, our approach 
yields a good approximation to the true distribution at the upper tail. Thus, the proposed series is more 
relevant to our problem than the series fitting method. 

Our approach to the upper tail approximation is based on the recent works by Raphaeli [22] and by 
Al-Naffouri and Hassibi [25]. First, let us explain Raphaeli's method. The procedure in Fig. 10 up to 
obtaining the MGF of the Gaussian quadratic form is common to both the sequence fitting method and 
Raphaeli's method. However, Raphaeli's method obtains the PDF by direct inverse Laplace transform of 
the MGF <&(s). Typically, the inverse Laplace transform of the MGF is represented as a complex contour 
integral and then the complex contour integral is computed as an infinite series by the residue theorem. 
However, to obtain the cumulative distribution function (CDF), which is actually necessary to compute 
the tail probability, Raphaeli's method requires one more step, the integration of the PDF, to obtain the 
CDF since the MGF <&(s) is the Laplace transform of the PDF. 

To obtain the CDF of a general Gaussian quadratic form, we did not use the MGF 3>(s), which is 
a bit complicated and requires an additional step, like Raphaeli, but instead we directly used a simple 
contour integral for the CDF (12), obtained by Al-Naffouri and Hassibi [25]. 6 Then, the contour integral 
was computed as an infinite series by the residue theorem. (Using the residue theorem is borrowed from 

5 In the case of the problem considered in [23], the outage defined in [23] is associated with the lower tail of the distribution 
and thus the series fitting method is well suited to that case. However, our system setup and considered problem are different 
from those in [23]. 

6 In [25], Al-Naffouri and Hassibi obtained the contour integral, (12) for the CDF of a Gaussian quadratic form. However, 
they did not obtain closed-form series expressions for the contour integral in general cases except a few simple cases. The main 
goal of [25] was to derive a nice and simple contour integral form for the CDF. 
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(a) (b) 

Fig. 11. Series fitting method versus direct inverse Laplace transform method: number of variables = 4, fj, = 0.51, Q = 
[1,0.5, 0,0; 0.5, 1,0,0; 0,0, 1,0; 0,0, 0,1], and £ = 0.31. (a) /3 = 1 and (b) = 2. (J) is the control parameter for the 
Laguerre polynomials in (44).) Note that the convergent speed of the series fitting method based on the Laguerre polynomials 
depends much on /3. In the case of /3 — 2, the series fitting method based on the Laguerre polynomials yields large errors at the 
upper tail. It is not simple how to choose f3 and an efficient method is not known. (One cannot run simulations for empirical 
distributions for all cases.) The series fitting method based on the power series shows bad performance, and it cannot be used 
in practice. 



Raphaeli's work.) Thus, our result is simpler than Raphaeli's approach and does not require the integration 
of a PDF for the CDF. 

As mentioned already, the series expansion in this paper has a particular advantage over the series 
fitting method considered in [23] for the outage event defined in this paper; The series in this paper fits 
the upper tail of the distribution well with a few number of terms. We shall provide a detailed proof for 
this in a special case in the next subsection. Thus, our series expressions for outage probability in MIMO 
interference channels are meaningful and relevant. 
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Appendix C 

Computational Issues and Convergence of the Obtained Series 
A. Computing higher order derivatives 

The general outage expression in Theorem 1 is given by 

Pr{outage} = Pr{log 2 (l + SINR^ m) ) < R™} 

» 



i=l 1 n=K ; — 1 

where 



PYn / V (s-l/Aj)A p |-v(9)|2 



ft ( s ) = ■ — TTx 7 • (47) 

s " 1/A * (l + 0» " 1/A,)A P ) 

To compute (46), we need to compute 

• {^i} (the eigenvalues of the Kd x ETd covariance matrix I] = SI/Avl/^), 

• {x| } ( tne elements of vector x = A" 1 / 2 *^/!, where /i is the mean vector of the Gaussian 
distribution), 

• and the higher order derivatives of gi(s). 

The computation of {A^} and {\\ } is simple since the sizes of the mean vector and the covariance 
matrix are Kd and Kd x Kd, respectively. Furthermore, the higher order derivatives of gi{s) can also 
be computed efficiently based on recursion [26], [22]. Note that 9i (s) = e los3i ^. Thus, the derivative of 
<7i(s) can be written as 

n-l , _\\ 

(s) = £ ( n 1 K } W P°g * (^] (n "° , n > 1 (48) 

where gf\s) and [log<7i(s)]0 denote the l-th derivatives of g-i(s) and log ft (s), respectively. Here, 
[log gi(s)]( n > can be computed from (47) as 

flog - r Sl i-m-ir- 1 T nlj-ir^X; - (n-l)!(-irSAg 

where <5i n is Kronecker delta function. Thus, for given gi{s) and [log ft (s)]^, we can compute gf\s) 
efficiently in a recursive way, as shown in (48). 
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B. Convergence analysis 

In this subsection, we provide some convergence analysis on the derived series expansion in Sec. III. 
Consider the general result in Theorem 1 for the CDF of a Gaussian quadratic form: 

1=1 1 n=Ki — \ \ / 

(49) 

where 

e»» exp i i^ v+i 1+(s _ 1/Xt)Xp 2^ q =i \X P I 

9i{s,y) 

Here, we explicitly use the variable y as an input parameter of the function gi{s) for later explanation. 
<?i"^( s >2/) denotes the n-th partial derivative of gi(s,y) with respect to s. (Here, k is the number of 
distinct eigenvalues of the Kd x Kd covariance matrix S and Kj is the geometric order of eigenvalue 
Aj. J2i=i K i = Kd.) The residual error caused by truncating the infinite series after the first N terms is 
given by 

RN(y) _ E ^_ E 2. 9f ) ( „, , (50) 

8=1 1 n=AT+l v \ / 

and we have 

Pr{y < y; infinite sum} = Pr{Y < y; truncation at N} + Rjy(y). 
The truncation error Rn{u) can be expressed as 

K 

R N (y) = Y J R i N(y), (51) 

i=i 

where 

for each 1 < i < k. Then, the magnitude of each term \R l N (y)\ in the truncation error is bounded as 

i /E,-=ilXi 



A » 1 VAi j=l /J n=JV+l 



(i),2\ n ~ K * +1 



(re - Ki + 1)! I Aj 



(53) 

As seen in Fig. 11, our series expansion fits the upper tail distribution first. Now, to assess the overall 
convergence speed of our series, for the same step as in Fig. 1 1, we ran some simulations to obtain an 
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empirical distribution, and computed the overall mean square error (MSE) between the truncated series 
and the empirical distribution over < y < 10 as 

1 200 2 

CDF MSE = ^2 Pr{^ < yf, N, type of series} - Pr{Y < yf, empirical} , 

200 I 

where {y{\ are the uniform samples of [0,10]. Fig. 12 shows the CDF MSE of the three methods in 
Fig. 11: the proposed series, the series fitting method with /3 = 1 and the series fitting method with 
(3 = 2. It is seen in Fig. 12 that the overall convergence of the proposed series can be worse than 
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Fig. 12. CDF MSE of the CDFs in Fig. 11 

the series fitting method at the small values for the number of summation terms for the setting in Fig. 
11. The bad overall convergence is due to worse fitting at the lower tail of the distribution, but the bad 
lower tail approximation is not important to our outage computation. (Please see Fig. 11.) Fig. 13 shows 
another case. In this case, the proposed series outperforms the series fitting method both in the overall 
convergence and in the upper tail convergence. It is seen numerically that the proposed series fits the 
upper tail distribution first. Now, we shall prove this property of the proposed series. However, it is a 
difficult problem to prove this property in general cases. Thus, in the next subsection, we provide a proof 
of this property when the number of distinct eigenvalues of the covariance matrix X is one, e.g., in the 
i.i.d. case. 

1) The identity covariance matrix case: Suppose that there is only one eigenvalue, A (> 0), with 
multiplicity k for the covariance matrix S. This case corresponds to Corollary 4, and the outage probability 
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Fig. 13. number of variables = 4, fl = 0.51, Q = I, and £ = [0.2641 0.0328 0.1963 0.1140; 0.0328 0.6097 - 
0.1739 0.1708; 0.1963 - 0.1739 0.8746 - 0.0022; 0.1140 0.1708 - 0.0022 0.1250]. In this case eigenvalues are 1.0000, 
0.6318, 0.2158, and 0.0259 with /3 = 1. (a) CDF, (b) CDF MSE. Uniform sample of y is taken over [0, 15.9], 



is given by 

n=K—l 

where 



pr { y < y} = i + ^^ eX p (--) ^ ^ ) (°^) re!(n - K + i)r ^ 



g{s,y) = z — w (55) 



and r/ 2 = J2j=i \ x^\ 2 - The residual error caused by truncating the infinite series after the first N terms 
is given by 

ri=7V+l v y 

Before we proceed, we first obtain the ra-th derivative of g(s, y) at s = 0, which is given in the following 
lemma. 

Lemma 1: For n > 0, 

n | 

5^(0, y) = -AV^-TXT A V^- (57) 
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Proof: Proof is given by induction. The validity of the claim for n = 0, 1 and 2 is shown by direction 
computation: 



5 (1) (o,y) 
<? (2) (o,y) 



ye 



s-l/X 



0! 



fc=0 



(0-fc) 



A fc y°- fc , 



ye^(s - 1/A) - e 



1! 



-A(y + A) = -A^^— ^Ay~ fc , 



s=0 fc=0 
(yeW(y a - y/X - 1) + e^y) (a - ±) 2 - 2e^(y S - y/X - l)(s - ±) 



(s - 1/A)* 

A(y 2 + 2Ay + 2A 2 ) = -A £ .^^aV^ 



s=0 



fc=0 



(2-fc)! 



Now, suppose that (57) holds up to the (n — l)-th derivative of g(s,y). From the recursive formula in 
(48), y (n) (0, y) is obtained as 

9 H (0,y) 



n-l / _ -i\ 

E ( k 5 {fc) (0,y)(log 5 (0, y ))("- fc ) 

fc=0 ^ ' 

n ^ 1 ^(°)(0,y)(log 5 (0,y))M + ( ra ~ 1 ^( 1 )(0,y)(log 5 (0,y))("- 1 ) 

+ (^lj)5 (n - 1) (o,y)(iogff(o,y)) (1) - 



+ ... 



(58) 



Since [logy(s)] = ys — log(s — 1/A), we can easily see that [log (7(0)] t 1 ) = y + X and [logy(0)]( n ) = 
(n — l)!A n for n > 2. Therefore, (58) can be rewritten as 

gW(0, y) =(n - l)\g(0,y)\ n + (n - l)g^(0,y)(n - 2)\X n ~ 1 +{^~ ^^(O, y)(n - 3)!A"" 2 + • • • 
+ (n - l)g^ (0, y) A 2 + y^" 1 ) (0, y)(y + A) 
=(n - 1)1(7(0, y)A n + (n - l)!y« (0, y^"" 1 + ^^V 2 > (0, y)A"" 2 + • • • 



+ (n - 1) 5 ("~ 2 ) (0, y) A 2 + A^" 1 ) (0, y) + yy^ 1 ) (0, y) 



(a) 



A 



.z=o 

n-l 



vfc=0 



(Z-fe)! 



E^ Ett-^av- 



m=0 
n-l 



(n — m — 1)! 



(Z-fe)! 



i _l_ ^ ^ ( ra jjj X m y n ~ m 



m=0 



(n — m — 1)! 



(59) 



_Z=0 \fc=0 

where (a) holds since (57) holds for all y(°)(0, y), ■ ■ ■ , g^ n ~^ (0, y) by the induction assumption. 



Here, consider the coefficient of each y 1 in (59) for i = 0, 



n. 
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i) y n is obtained only when m = 0. The coefficient of y n from (59) is therefore given by —A. It 
corresponds to the coefficient of y n in (57). 

ii) For < p < n, the coefficient of y n ~P is obtained by considering all (/, k) that satisfies I — k = 
n — p due to the first term in the right-hand side (RHS) of (59), and m = p due to the second term 
of the RHS of (59). In the first case, we obtain y n ~P with the following pairs (I, k) = (n — l,p— 1), 
(n — 2,p — 2), • • • , (n — p, 0). For these (I, k) pairs, we have 



-a E = E 



(60) 



In the second case of m = p, we have 



in - l)' 

TwW _p . (61) 

[n — p — iy. 



Finally, the coefficient of y n 9 is given by adding (60) and (61): 

-x( (n " 1)! +P i^Vv-p 

l (n — p — 1)! (n — p)\ J 

= - \ (w " 1) lJ i+ — 

(n — p — 1J! \ n — p I 

71 ' 



(n — p)! 

which is equivalent to the coefficient for y n_ P in (57) (0 < p < n). Thus, (57) holds for g^(0,y). 

■ 

Note that # (n) (0,y) < for all n > from (57). Therefore, i?jv(y) < for all N and y and |c/ (n) (0,y)| = 
-9 {n) (0,y). 

Now, consider the residual error term Rn{u) in (56). The magnitude of the residual error can be upper 
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bounded as follows: 

n=N+l V ' 

exp(-, ? 2 ) , n ~ „ (rj 2 / X) n ~ K+1 



• ex p(-l) E (-5 (n) (o,y)) : 



A* A^/ * ^""n!(n -«+l)! 



A- A^/ — + 1)1 

A K 6XP l AJ ^ n! 9 lU ' yj l2AJ (ra-K + 1)! 



n=JV+l v y 

n=JV+l 

i/c-1 , 00 1 ' 1 \n 



^ r exptf).ex P (-f) £ >>(0,,)(i 



n=7Y+l 

( 6 ) 2 K - X . o. / V\ ^ 1 . / 1 \ n 



A 

< ^-exp(r/ 2 )-exp(-|) " £ ^ <? (n) (0, 2/) 

( C ) / y 



exp(r/ ) ■ exp (--] g — ,y 



A ^ w y r V A7 3 V 2A 

(<*) 2 K - X 2 / exp(y/2A) 

■ exp(r/ j • exp I — — 



A ^' ' *\ XJ -1/2A 



= 2 K exp(r ? 2 ).exp(-^j (62) 

where (a) is from ^ < exp(7) = YlpLo l p /p ] - for an y 7 > 0, (b) is from the fact that summand is 
negative, (c) is by using the Taylor series expansion, and (d) is from (55). Since r\ is a fixed constant, 
from (62), for any N > 

lim |i?iv(2/)| = 0. (63) 

Thus, it is clear that the proposed series converges from the upper tail distribution! 

Now, let us consider the residual error magnitude as a function of y for given N. From (57), we have 

dgW(0,y) 



dy 



= ng^ n ~ 1 \o,y). (64) 
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Differentiating Rn(v) with respect to y yields 



n\(n — k + 1)! 

n=N+l v ; 



+ exp(-V) exp /_yx g 5 9 W(0,y) (r/ 2 /A)"~ K+1 



A 

Furthermore, from (57) we have 



X K V A/ ^— ' dy nUn — k + 1)1 

n=N+l y v ; 

n=N+l \ > ' \ / 



lgW(0, y )+ng( n - 1 \0,y) = y r 



(66) 



By substituting (66) into (65), we have 



dR N {y) exp(-?7 2 ) / y\ ~ (n 2 /X) n ~ K+1 y n 



<% A K 6 ^ V A^ ^— ' n!(n — k + 1)! 

tf n=iV+l V ' 

which is positive. Since Rn(v) < 0, lim y _ 5 . 00 Rn(v) = and dR Qy V ^ > 0, the residual error magnitude 
monotonically decreases as y increases and the maximum error occurs at y = for any given N. 
Now, let us compute the worst truncation error R^(0), which is given by 

A' 1 nun — k + 1 ! 

n=7V+l V y 

From (57), we have g^(0,0) = — n!A n+1 . Therefore, 



vi i*) y (fn . 



v ; X K ^ v ; n (n-K + l) 

n=iV+l v y 

eX p(- r? 2) - (?7 2 /A) n-« + l 



n=N+l V ' 

A K ^ + 

n=AT+l v y 

exp(-r? 2 ) £ ^J— - (69) 



From (54), N > k — 2. For general A r > k — 2, let m = n — k + 1. Then, 



R N (0) = -exp(- ?? 2 ) ^ 

TO=-/V-K+2 



??? 



Note that X]m=Af-K+2 ml * s tne res idual error of the Taylor series expansion of exp(x) after the first 
(JV — k + 1) terms. By the Taylor theorem, 

£ ^= i" } ^, ex P (a» 2 ) (70) 

m=JV-«;+2 v y 
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where some a G [0, 1]. Therefore, the worst truncation error is given by 

. . ( n 2\N-K+2 ( n 2\N-tt+2 

|«v(0)| = exp ((„ - V) x ( ^ K + 2) , < ( ^, + 2)| . PI) 

where the inequality holds since exp((a — l)^ 2 ) < 1 for < a < 1. Furthermore, the residual error 
magnitude is a strictly decreasing function of N for any y, 



\R N (y)\ > \Rn+i(v)\. (72) 



This can be shown easily as follows. 



exp(-r ? 2 ) / ^ Mfn , (rj 2 / X) n ~ K+1 



vy; A K V A/ ' v ,yj n\(n - k + 1)! 

_ exp(-, 2 ) / yx f ^ (n) (.VA)-^ (jV+1) (^/A)^ ) 

=Rn+i{v) + ^^ exp r aJ • 5 (0 ^ ) (jv +1 ) !(i v_ K+2)! - 

Since R N {y) <0 and y (7V+1) (y) < for all y > and iV, we have (72). Now, based on (71) and (72), 
with given x fc an ^ °h> we can compute the required number iV of terms in the series to achieve the 
desired level of accuracy since rj 2 is known. 




Finally, consider the worst case of N = k — 2 and y = 0: 

i^ 2 (0) = -exp(- ?? 2 ) £ (r] ) =-eM-if)Y. - lU r = - 1 

z — ' (n — K + l)! m! 

n=K—l m=0 
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where the second equality is by replacing m = n — k + 1. It is easy to see that the worst case error is -1 
in the identity covariance matrix case. Fig. 14 shows the performance of the proposed series expansion 
in the case of the identity covariance matrix. The numerical results well match our theoretical analysis 
in this subsection. From the figure, it seems reasonable to choose ./V > 20 ~ 30 for accurate outage 
probability computation. 
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